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ABSTRACT 

A number of studies have shown that the convective stability criterion for the in- 
tracluster medium (ICM) is very different from the Schwarzchild criterion due to the 
effects of anisotropic thermal conduction and cosmic rays. Building on these studies, 
we develop a model of the ICM in which a central active galactic nucleus (AGN) ac- 
cretes hot intracluster plasma at the Bondi rate and produces cosmic rays that cause 
the ICM to become convectively unstable. The resulting convection heats the intra- 
cluster plasma and regulates its temperature and density profiles. By adjusting a single 
parameter in the model (the size of the cosmic-ray acceleration region), we are able to 
achieve a good match to the observed density and temperature profiles in a sample of 
eight clusters. Our results suggest that convection is an important process in cluster 
cores. An interesting feature of our solutions is that the cooling rate is more sharply 
peaked about the cluster center than is the convective heating rate. As a result, in sev- 
eral of the clusters in our sample, a compact cooling flow arises in the central region 
with a size r^f that is typically a few kpc. The cooling flow matches onto a Bondi flow 
at smaller radii. The mass accretion rate in the Bondi flow is equal to, and controlled 
by, the rate at which mass flows in through the cooling flow. Our solutions suggest that 
the AGN regulates the mass accretion rate in these clusters by controlling r^: if the 
AGN power rises above the equilibrium level, r^f decreases, the mass accretion rate 
drops, and the AGN power drops back down to the equilibrium level. 

Subject headings: cooling flows — galaxies:clusters:general — galaxies: active — 
convection — magnetic fields — turbulence 
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1. Introduction 

Active galactic nuclei (AGNs) have enormous mechanical and radiative luminosities. If an 
AGN's power can be transferred to the surrounding interstellar and intergalactic media, the result- 
ing heating can have a large effect on the ambient plasma. There has recently been great interest 
in this process of "AGN feedback," its role in galaxy formation, and the possibility that AGN feed- 
back solves the over-cooling problem (Suginohara & Ostriker 1998, Lewis et al 2000, Tomatore 
et al 2003, Nagai & Kravtsov 2004) and cooling-flow problem for clusters of galaxies (Bohringer 
et al 20001; David et al 2001; Tamura et al 2001; Molendi & Pizzolato 2001; Blanton, Sarazin, & 
McNamara 2003; Peterson et al 2001, 2003). 

One of the main unsolved problems for AGN feedback is to understand how AGN power is 
transferred to the diffuse ambient plasma. A number of mechanisms have been investigated, in- 
cluding Compton heating (Binney & Tabor 1995; Ciotti & Ostriker 1997, 2001; Ciotti, Ostriker, 
& Pellegrini 2004, Sazonov et al 2005), shocks (Tabor & Binney 1993, Binney & Tabor 1995), 
magnetohydrodynamic (MHD) wave-mediated plasma heating by cosmic rays (Bohringer & Mor- 
fiU 1988; Rosner & Tucker 1989; Loewenstein, Zweibel, & Begelman 1991), and cosmic-ray 
bubbles produced by the central AGN (Churazov et al 2001, 2002; Reynolds 2002; Bruggen 2003; 
Reynolds et al 2005), which can heat intracluster plasma by generating turbulence (Loewenstein 
& Fabian 1990, Churazov et al 2004) and sound waves (Fabian et al 2003; Ruszkowski, Bruggen, 
& Begelman 2004a,b) and by doing pdV work (Begelman 2001, 2002; Ruszkowski & Begel- 
man 2002; Hoeft & Bruggen 2004). Despite this substantial progress, it is still not clear how AGN 
feedback controls the density and temperature profiles of the ambient plasma in a way that is both 
self -regulating and consistent with observations. 

In this paper, we focus on clusters of galaxies and explore the hypothesis that central AGNs 
heat and regulate the intracluster plasma by causing the intracluster medium to become convec- 
tively unstable, a scenario that was investigated in two earlier studies [Chandran (2004) (hereafter 
Paper I) and Chandran (2005) (hereafter Paper II)]. At first glance, this hypothesis seems obviously 
incorrect, since observations show that the specific entropy s in intracluster plasmas increases with 
radius r. However, several recent studies have shown that the Schwarzchild criterion {dsjdr > 0) 
does not apply to low-density, magnetized plasmas such as those found in clusters, in which the 
charged-particle gyroradii are much less than the Coulomb mean free path. In such plasmas, heat 
and charged particles diffuse primarily along magnetic field lines, and only weakly across the mag- 
netic field. This anisotropy turns out to have a profound effect on convective stability, as shown 
analytically by Balbus (2000, 2001) and numerically by Parrish & Stone (2005, 2007). These au- 
thors considered a stratified plasma in which the gravitational acceleration is in the —z direction and 
the equilibrium magnetic field is in the xy-plane and showed that the convective stability criterion 
is dT/dz> 0, not ds/dz> 0, where T is the temperature. When cosmic rays are present, the con- 
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vective stability criterion becomes nkBclT /dz + dpcr/dz > 0, as shown analytically by Chandran 
& Dennis (2006) and numerically by Rasera & Chandran (2007). Here, n and pcr are the thermal- 
plasma number density and cosmic-ray pressure, respectively. In galaxy clusters, the gravitational 
acceleration is in the —r direction, and the convective stability criterion is 



dT dpci 

nkB— + ——>0. (1) 
dr dr 



(Paper II and appendix |B] provide a more extensive discussion.) Although dT /dr > in cluster 
cores, equation ([T]) shows that cosmic rays produced by an AGN at the center of a cluster can lead 
to convective instability, since centrally produced cosmic rays satisfy dpc^/dr < 0. 

In this paper, we construct a spherically symmetric, steady-state model of convective intra- 
cluster plasmas using mixing-length theory, and compare this model to observations. We assume 
that a central supermassive black hole accretes hot intracluster plasma at the Bondi rate, and con- 
verts a small fraction of the accreted rest-mass energy into cosmic rays that are accelerated by 
shocks within some distance rgource of the center of the cluster. The resulting cosmic-ray pressure 
gradient leads to convection, which in turn heats the thermal plasma in the cluster core by advecting 
internal energy inwards and allowing the cosmic rays to do pdV work on the thermal plasma. The 
model also includes thermal conduction, cosmic-ray diffusion, and radiative cooling. The model 
involves much less emission from plasma at temperatures below one-third of the cluster's average 
temperature than the cooling flow model (Fabian 1994), and thus offers a possible solution to the 
cooling-flow problem. 

We compare the density and temperature profiles predicted by the model to the profiles in- 
ferred from X-ray observations of eight clusters. We adjust a single parameter, the size rgource of 
the cosmic-ray acceleration region, to optimize the fit. The model solutions match the observa- 
tions well, with the exception of the density with the central ~ 50 kpc of Sersic 159-03, which is 
underestimated by the model. We suggest a possible explanation for this discrepancy in section[3l 
We also find that the cosmic-ray luminosities of the AGN in our sample are strongly correlated 
with the observationally inferred mechanical luminosities of these AGN. Our results suggest that 
AGN-driven convection is an important process in cluster cores. 

An attractive feature of this model and other models based on AGN feedback and Bondi 
accretion is that they are self-regulating. One argument for why Bondi accretion is self -regulating 
was advanced by Nulsen (2004) and Bohringer et al (2004), who noted that the Bondi accretion 
rate is a monotonically decreasing function of the specific entropy near the center of the cluster. 
Thus, if the central plasma becomes too cool, the Bondi accretion rate rises, the AGN feedback 
heating increases, and the specific entropy of the central plasma rises back to its equilibrium value. 
In this paper, we offer an additional explanation for how AGN heating on large scales (> 5 kpc) 
can regulate the mass accretion rate onto the central black hole. In our solutions, we find that the 
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radiative cooling rate is more sharply peaked about the center of a cluster than is the convective 
heating rate. As a result, in several of the clusters in our sample, the central region becomes a 
cooling flow. The radius of this cooling flow, r^f, is typically a few kpc in our solutions. At smaller 
radii, the flow makes a transition from a cooling flow to a Bondi flow. However, as in the work 
of Quataert & Narayan (2000), the mass accretion rate of the inner Bondi flow is controlled by 
the surrounding cooling flow. In our model, which has no mass dropout, the mass accretion rate 
is approximately the plasma mass interior to r^f divided by the cooling time at r^f. The AGN then 
regulates the mass accretion rate by controlling rcf : if the AGN power rises above the equilibrium 
level, the size of the central cooling flow decreases, the mass accretion rate drops, and the AGN 
power then drops back down to the equilibrium level. 

This paper extends the previous models of paper I and paper II in several ways. In contrast to 
paper I, the present paper takes into account the role of anisotropic thermal conduction and cosmic- 
ray diffusion, which strongly modify the convective stability criterion. In contrast to paper II, we 
take the cosmic -ray acceleration to occur within a relatively small fraction of the total volume at 
any given radius, which allows for localized pockets of excess cosmic -ray pressure that tend to rise 
buoyantly. We also take into account the nonzero average radial velocity, and compare the model 
to a larger sample of clusters. 

The rest of this paper is organized as follows. We present the basic equations of the model 
in section |2l In section [3l we compare our model calculations to observations. In section |4] we 
consider the radial profiles of the different heating rates and the factors that determine whether 
AGN feedback or thermal conduction is the dominant heat source at r < 100 kpc. In section[5]we 
discuss the central cooling flows that arise in our model solutions for several of the clusters in our 
sample. We also comment in section [5] on the viability of the Bondi accretion model for the AGN 
at the centers of clusters. We summarize our results in section |6l We present results on the radial 
profiles of the turbulent velocity and cosmic -ray pressure in appendix|Al In appendix|B]we present 
a systematic derivation of the two-fluid mixing-length theory that we employ in our model. 



We describe the intracluster medium using a standard set of two-fluid equations for cosmic 
rays and thermal plasma (Drury & Volk 1981, Jones & Kang 1990), modified to include thermal 
conduction, viscous dissipation, and radiative cooling: 



2. Model equations 



(2) 



dv 



v(p+Pcr)-pv$-v-n 



(3) 




vise 5 
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^ = -ypV ■ V + (Y- 1 ) [//vise + V ■ (K ■ Vr) - /?] , (4) 

and 

^ = -Ycri^crV ■ V + V ■ (D ■ Vpcr) + (Ycr " l)^source, (5) 

where 

p is the plasma density, v is the bulk velocity of the two-fluid mixture, p and pcr are the plasma 
and cosmic-ray pressures, T is the plasma temperature, 4> is the gravitational potential, Ilvisc is the 
viscous stress tensor, y and Ycr are the plasma and cosmic-ray adiabatic indices (which are treated 
as constants), //vise is the rate of viscous heating, K is the thermal conductivity tensor, R is the 
radiative cooling rate, ^source is the rate of injection of cosmic -ray energy per unit volume by the 
central radio source, and D is an effective momentum-averaged cosmic -ray diffusion tensor. For 
the calculations presented in section [3l we set y= 5/3 and Ycr = 4/3. We ignore radiative cooling 
of cosmic rays, which is reasonable if protons make the dominant contribution to the cosmic -ray 
pressure. We also neglect Coulomb interactions between cosmic rays and thermal plasma, as well 
as wave-mediated heating of the thermal plasma by cosmic rays. As discussed below, we take 
into account the effects of the magnetic field on K and D, but we neglect the Lorentz force and 
resistive dissipation. In the following subsections, we describe the approximations we use to solve 
the above equations. 



2.1. Mixing length theory 

To account for convection, we write each fluid quantity as an average value plus a turbulent 
fluctuation: 

v=(v) + 6v, (7) 

p = {p)+hp, (8) 

etc, where (...) denotes an average over the turbulent fluctuations. We take the averaged quantities 
to be spherically symmetric and independent of time, and we treat the fluctuating quantities as 
small. To obtain equations for the average cluster properties, we average equations Q through ([5]). 
We evaluate the averages (pv), (yp), and (vpcr) in equations Q, ©, and ^ using a two-fluid 
mixing-length theory that we describe in appendix |Bl The essential idea behind this theory is that 
the amplitudes of the turbulent fluctuations increase as the average plasma and cosmic -ray profiles 
move past the point of marginal stability towards increasing degrees of convective instability. As 
this happens, the magnitudes of the internal energy flux (vp) / (y— 1 ) and the cosmic-ray energy flux 
{vpcx) I (Ycr — 1) increase, which in turn affects the density, temperature, and cosmic-ray-pressure 
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pro files. The two-fluid mixing length theory provides an approximate way of determining the 
resulting profiles as well as the r-dependent turbulent velocity in a self-consistent way. A key 
parameter of the model is the mixing length /, which characterizes the length scale of the convective 
turbulence. We set 

I = OAr. (9) 



2.2. Hydrostatic equilibrium 

We assume that the convection is subsonic and confine our model to r > 0.2 kpc, so that 
the average radial velocity remains subsonic throughout our solutions. As a result, we can to a 
reasonable approximation drop the inertial terms in the average of equation (|3]). The viscous term 
in equation (|3]) is important primarily for dissipating small-scale velocity fluctuations and can also 
be neglected in the average of equation ([3]). The average of equation (|3]) then reduces to 

^(Ao.) = (10) 

where 

Ptot=P+Pci- (11) 



2.3. Gravitational potential 

We take the gravitational potential to be the sum of four components, 

4> = 4>c-f4>s+$bh + ^p, (12) 

where 4>c is the contribution from the the cluster's dark matter, 4>s is the contribution from the stars 
in the brightest cluster galaxy (BCG), 4>bh is the contribution from the black hole at r = 0, and 4>p 
is the contribution from the intracluster plasma. We take the cluster dark matter to have an NFW 
density profile (Navarro, Frenk, & White 1997), 

PDM = — ] 79-5 (13) 

where 

_ 200 

'~ 3 [ln(l + c) - c/(l + c)]' ^ ^ 

rs is the scale radius, c is the concentration parameter, and Pcdt = SH^/SnG is the critical density 
at the redshift z of the cluster. The latter is calculated assuming Q.o = 0.3, Q.^ = 0-7, and Hq = 
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Table 1 : Parameters used in determining the gravitational potential 



Cluster 


BCG 


rs 


c 


z 


Mb 


B-V 


Lb 


Re 


Mbh 






(kpc) 










(1O"Lb,0) 


(kpc) 


(lOX?) 


Virgo 


NGC 4486 (M87) 


560 


2.8 


(see below) 


-21.96 


0.93 


0.938 


5.03 


1.38 


Abell 262 


NGC 0708 


85 


8.62 


0.0155 


-21.08 


1.06 


0.417 


25.6 


0.555 


Sersic 159-03 


ESQ 291-009 


159 


6.56 


0.0572 


-22.16 


1.00 


1.13 


29.5 


1.92 


Abell 4059 


ESQ 349-010 


744 


2.7 


0.0466 


-22.73 


1.06 


1.91 


24.5 


4.12 


Hydra A 


PGC 026269 


77 


12.3 


0.0550 


-22.97 


0.82 


2.38 


39.6 


4.12 


Abell 496 


PGC 015524 


129 


7.75 


0.0322 


-22.48 


1.12 


1.51 


49.9 


3.27 


Abell 1795 


PGC 049005 


430 


4.21 


0.0639 


-22.04 


1.00 


1.01 


40.3 


1.66 


Perseus 


NGC 1275 


481 


4.09 


0.0179 


-22.62 


0.53 


1.72 


15.3 


1.89 



The NFW parameters and c describe the clusters' dark matter density profiles. For Virgo and c are taken from 
McLaughlin (1999). For Hydra A, r, and c are taken from David et al (2001). For all other clusters, and c are taken 
from table 1 of Piffaretti et al (2005). Redshifts z are taken from Kaastra et al (2004), except for Virgo — Kaastra 
et al (2004) take the distance to Virgo to be 16 Mpc, and we use the same value. Absolute B-band magnitudes Mb 
and B — V color indices for the brightest cluster galaxies (BCGs) are taken from the "Hyperleda" database of Paturel 
et al (2003). The BCG effective radii Re are taken from Schombert (1987) for Perseus and Abell 1795, from Graham 
et al (1996) for Hydra A, Abell 262, and Abell 496, and from "Hyperleda" for Virgo, Sersic 159-03 and Abell 4059. 
Lb is the BCG B-band luminosity. The black hole masses are determined using the mass-luminosity relation given in 
equation (6) of Lauer et al (2007). 

70 km s^^Mpc^^ The values of r^, c, and z for the eight clusters we consider in section[3]are taken 
from the literature and listed in table [TJ 

We take the stellar mass density to have a Hernquist profile in which the stellar mass interior 
to radius r is 

Mstars r = (15) 

where Mq is the total stellar mass and a is a scale length equal to i?e/1.8153, where Rg is the 
radius of the isophote enclosing half the galaxy's light. (Hernquist 1990) As in Graham et al 
(2006), we set Mq = TbLb, where Lb is the BCG B-band luminosity, and Tg = 5.3Mq/Lb,q is the 
B-band stellar mass-to-light ratio for a 12-Gyr-old single stellar population (Worthey 1994). We 
set L5/L5 = lO^-'^i^B/^j-MB) ^ where Mb.q and Mb are, respectively, the solar and BCG absolute 
B-band magnitudes, and M5 © = 5.47 (Cox 2000). The values of Rg and Mb for each cluster are 
taken from the literature (see table [I]). 
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We determine the black-hole mass using the mass-luminosity relation given in equation (6) of 
Lauer et al (2007): 

logf^ ) = 8.67-0.528(Mk + 22), (16) 

where My is the BCG absolute V-band magnitude. We set My = Mb — (B — V), where Mb and the 
B — V color index for each cluster are taken from the "Hyperleda" database (Paturel et al 2003) 
and listed in table [TJ The resulting values of Mbh for each cluster are also listed in table [IJ 

The contribution to the gravitational potential from the intracluster plasma 4>p is not deter- 
mined ahead of time, but is instead obtained by solving V^^p = — 47cG(p) , where (p) is the average 
plasma density that results from solving the model equations. 



2.4. Radiative cooling and chemical composition 

We use the analytic fit of Tozzi & Norman (2001) to approximate the full cooling function for 
free-free and line emission: 



R = mrie 



0.0086 (-^^ +0.058 f-^^ +0.063 
1 keV / V 1 keV 



• 10^22 ergs cm^ s"\ (17) 



where is the ion density, is the electron density, is the Boltzmann constant, and the numer- 
ical constants correspond to 30% solar metallicity. Because we treat the turbulent fluctuations as 
small, we can replace rie, n-i, and T in equation (flTI) by their average values when calculating {R). 
We take the intracluster plasma to be fully ionized and to have a uniform chemical composition, 
with a hydrogen mass fraction of X = 0.7 and a helium mass fraction Y = 0.29. We take the metals 
to have a mean charge to mass ratio equal to that of helium. The mean molecular weight is then 

fu=- = 0.62. (18) 

(ne + HiJm// 

The mean molecular weight per electron is then 

A^e = -^ = 1.18. (19) 

In addition, 

— = 0.91, (20) 

and 

rip 

— = 1.21, (21) 
nn 

where nn is the hydrogen number density. 
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2.5. Transport 

Cluster magnetic fields are easily strong enough to cause cosmic rays and heat to diffuse 
primarily along magnetic field lines, so that 

K~K||^, (22) 

and 

D~Z)||^, (23) 

where b is the magnetic field unit vector, and K|| and D|| are the parallel conductivity and diffusivity. 
We take the parallel conductivity to be the classical Spitzer thermal conductivity (Spitzer & Harm 
1953,Braginskii 1965), 

where In Ac is the Coulomb logarithm. The local anisotropy of K and D turns out to be critical 
for convective stability, as discussed by Balbus (2000,2001), Parrish & Stone (2005,2007), Chan- 
dran & Dennis (2006), and Rasera & Chandran (2007), and we take this anisotropy into account 
in our mixing length theory for intracluster convection. [See, e.g., the discussion preceding equa- 
tion (|B29I) .1 However, when we average equations ^ through ([5]) and solve for the structure of 
the ICM, we are interested in the transport of heat and cosmic rays over distances much greater 
than the correlation length of the magnetic field, Ig, which is ~ 1 — 10 kpc (Kronberg 1994; Taylor 
et al 2001, 2002; Vogt & Ensslin 2003, 2005 - see Schekochihin et al 2006 and Schekochihin & 
Cowley 2006 for a recent discussion of intracluster magnetic fields and turbulence). For transport 
over such large scales, averaging over the turbulent magnetic field leads to an effectively isotropic 
conductivity, which we denote Kj, that is reduced relative to Ky (Rechester & Rosenbluth 1978, 
Chandran & Cowley 1998). Theoretical studies find that the reduction is by a factor of ~ 5 — 10 
(Narayan & Medvedev 2001, Chandran & Maron 2004, Maron, Chandran, & Blackman 2004). In 
this paper, we assume that 

Kr = ^. (25) 



(v.(K-vr))- 2 



1 d 
dr 



2, d 



dr 



We take the average of the conductive heating term to be given by 

(26) 

with T set equal to {T) in equation (l24l) . Similarly, we assume that 

(27) 



(V-(D-Vper))- 2 



1 d 
r^ dr 



r Dci—{pci) 
dr 
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We take the value of Drv to be 



(28) 



where Do = 10^^ cm^/s and Vd = 10 km/s. The vj term is loosely motivated by a simplified pic- 
ture of cosmic-ray "self-confinement," in which cosmic rays are scattered by waves generated by 
the streaming of cosmic rays along field lines. If, contrary to fact, the field lines were purely ra- 
dial, efficient self-confinement would limit the average radial velocity of the cosmic rays to the 
Alfven speed va, allowing the cosmic rays to travel a distance r in a time ~ r/vA- For constant 
va, this scaling can be approximately recovered by taking the cosmic rays to diffuse isotropically 
with Dci <^ r, the scaling that arises from equation (l28l) when v^r ^ Dq. This self-confinement 
scenario is too simplistic, since in clusters field lines are tangled, va varies in space, and it is not 
known whether cosmic rays are primarily scattered by cosmic -ray-generated waves or by magneto- 
hydrodynamic (MHD) turbulence excited by large-scale stirring of the intracluster plasma. It is not 
clear, however, how to improve upon equation (|28|) . Self-confinement in the presence of tangled 
field lines is not well understood, and the standard theoretical treatment of scattering by MHD tur- 
bulence, which takes the fluctuations to have wave vectors directed along the background magnetic 
field, is known to be inaccurate (Bieber et al 1994, Chandran 2000, Yan & Lazarian 2004). A more 
definitive treatment must thus await further progress in our understanding of MHD turbulence and 
cosmic-ray transport. The value of is needed in the mixing length theory developed below. We 
assume that Da/D^^ = Kx/^p and thus set 



2.6. The mass accretion rate and cosmic-ray luminosity of the central AGN 

We assume that the black hole at r = in our model, with a mass Mbh given by equation (fT6l) , 
accretes intracluster plasma at the Bondi (1952) rate. 



where Cs is the adiabatic sound speed, and p and Cs are evaluated using the average plasma pa- 
rameters at the radius ri = 0.2 kpc, which defines the inner boundary of our model solutions. We 
assume that this accretion powers a jet that leads to shocks, which in turn accelerate cosmic rays. 
We take the cosmic -ray luminosity to be 



D I = 8A 



(29) 



M 




(30) 



'cr — 



(31) 



where 



ri = 5x 10~3. 



(32) 



- 11 - 



An argument against Bondi accretion in clusters is that the radiative luminosities of AGNs in 
elliptical galaxies are typically several orders of magnitude smaller than the nominal Bondi accre- 
tion power, given by /'Bondi — 0- l^Bondi^^, where MBondi is the Bondi accretion rate given in equa- 
tion (1301) . (Allen et al 2006) However, the mechanical luminosities /^rnech 

of these AGN are often 

much larger than their radiative luminosities. Moreover, in a recent study of nine AGNs in nearby 
x-ray luminous elliptical galaxies, Allen et al (2006) found a strong correlation between PBondi (as 
calculated from the observed plasma temperature and density profiles) and L^ech (as inferred from 
the energies and time scales required to inflate the observed x-ray cavities). Allen et al (2006) 
found that Lmech can be related to /feondi by a power-law fit of the form log(PBondi/ 10^^ erg s^^) = 
ci +C2log(Lniech/10'^-^ erg s^i), with ci = 0.65 ±0.16 and C2 = 0.77 ±0.20, and that the fraction 
of Mc^ that is converted into mechanical luminosity ranges from 1.3% for a jet power of 10"*^ erg/s 
to 3.7% for a jet power of 10"*^ erg/s. Results consistent with these were also found by Tan & 
Blackman (2005). These authors reviewed studies of M87 and estimated that Lmech is about an or- 
der of magnitude larger than the radiative luminosity, and that Lmech ~ O.OlMBondiC^. Our choice 
of T] = 0.005 is smaller than the accretion efficiencies found in these studies, in part to provide a 
more conservative estimate, and in part because only part of the mechanical energy is converted 
into cosmic rays. 

We note that Bondi accretion in clusters has been considered previously by a number of au- 
thors (e.g., Quataert & Narayan 2000, Di Matteo et al 2002, Nulsen 2004, Bohringer et al 2004, 
Springel et al 2005, Cattaneo & Teyssier 2007). Also, in Tan & Blackman's (2005) analysis, part 
of the reason for the small value of r\ is that part of the mass flowing in through the Bondi radius 
never reaches the central black hole because it forms stars in a gravitationally unstable disk. Thus, 
the Bondi accretion rate in our model may be significantly higher than the time derivative of the 
mass of the central black hole. 

Pizzolato & Soker (2005) and Soker (2006) considered a different "cold feedback" scenario 
for mass accretion, in which cold gas fuels the central AGN. In section [5] we address several issues 
related to the question of whether one expects Bondi accretion or some form of cold feedback in 
clusters. 



2.7. Cosmic-ray acceleration by the central radio source 

The spatial distribution of cosmic-ray injection into the ICM is not precisely known. Some 
clues are provided by radio observations, which show that cluster-center radio sources (CCRS) 
differ morphologically from radio sources in other environments. As discussed by Eilek (2004), 
roughly half of the CCRS in a sample of 250 sources studied by Owen & Ledlow (1997) are "amor- 
phous," or quasi-isotropic, presumably due to jet disruption by the comparatively high-pressure. 
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high-density cluster-core plasma. With the exception of Hydra A, the CCRS in the Owen-Ledlow 
(1997) study are smaller than non-cluster-center sources, with most extending less than 50 kpc 
from the center of the host cluster (Eilek 2004). Given these findings, we take the cosmic -ray 
acceleration to be concentrated within the cluster core. 

In paper II, it was assumed that the cosmic rays are accelerated in an approximately volume- 
filling manner. In contrast, in this paper, it is assumed that cosmic -ray energy is injected into the 
intracluster medium in only a fraction of the volume at any given radius. We then take 

■^'source ~ (^source) "l~ S^source 5 (33) 

where 

(£source)=5oe-'''/''— (34) 

can be thought of as an average of ^source over spherical polar angles. The constant rsource is a 
free parameter that characterizes the size of the cosmic-ray acceleration region. The constant 5o is 
determined on energy grounds from the equation = 4-n Jq dr {E^omceif)) and equation (|3T| ). 
After determining (i'sourcel'"))^ we set 

S£'rms ~ ''12 (■^'source) 7 (35) 

where S^rms is the rms value of S^source^ and r\2 is a constant that is related to the volume filling 
factor of the cosmic-ray acceleration region. For example, suppose that ^source = C = constant 
in a fraction fa of the volume between radius r and r -|- dr, and that ^source = in the remainder 
of the volume between r and r + dr. In this case, (£'source('')) = fcrC, ([£'source]^) = /crC^, and 

5£rms = a/ {[E^ourceir) - {E^ourctir))]'^) = {Esomce{r)) \J fcr^ - 1. For the Calculations presented 
below, we set r\2 = 2.5, which corresponds to = 0.138. These fluctuations in the cosmic -ray 
source term drive fluctuations in the fluid quantities and contribute to convection. This effect is 
incorporated into the two-fluid mixing length theory presented in appendix |Bl The fluctuations in 
£^source result in larger fluctuations (spatial variations) in and p than in the model of paper II, 
which in some sense represent the "cosmic-ray bubbles" or X-ray cavities seen in about one-fourth 
of the clusters in the Chandra archive (Birzan et al 2004). 



2.8. Summary and numerical method 

The approximations described above lead to a set of coupled ordinary differential equations 
for the average density, temperature, and cosmic-ray pressure and the rms turbulent velocity. These 
equations are presented in appendix |Bl We solve this set of equations using a shooting method, in 
which we guess the electron density, temperature, and cosmic-ray pressure at the inner radius of 
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our model (ri =0.2 kpc) and then update these guesses until the model solution satisfies the three 
boundary conditions at the outer radius router- These outer boundary conditions are the observed 
electron density 

'^e, outer ^uid temperature Touter router? tiiid a condition on (c?/?cr/ dv^ at router? which 
amounts to requiring that (;?cr) ^ as r ^ oo. The value of router for a cluster is taken to be the 
radius of the first observational data point outside the cluster's cooling radius, rcooh given in tabled 
(except for Virgo, for which we take router to be the outermost data point, which lies inside of rcooi-) 
The values of router? "outer? and Touter are listed in tabled After finding the values of Ue, T , and pcx 
at r\ needed to match the boundary conditions at router? we integrate the equations out to radii 
greater than router as needed to compare to the data. A more extensive discussion of our numerical 
method is given in appendix |Bl 



3. Comparison to observations 

We compare our model solutions with observations of the central regions (r < 0.25rvir? where 
rvir = CTs is the virial radius) of eight clusters: Virgo, Abell 262, Sersic 159-03, Abell 4059, 
Hydra A, Abell 496, Abell 1795, and Perseus. Temperature and hydrogen-number-density {hh) 
profiles for these clusters are taken from table 5 of Kaastra et al (2004). Redshifts (z) and angular- 
diameter distances Jscdm are given in table 1 of Kaastra et al (2004). The data of Kaastra et 
al (2004) are obtained assuming a standard cold dark matter (SCDM) cosmology with f2 = 1 and 
Hq — 50 km s^^ Mpc^^. We convert to a ACDM cosmology with Q.q = 0.3, Q.a.o = 0.7, and 
Hq — 70 km s^^ Mpc^^ by calculating the ratio of angular-diameter distance in the two cosmolo- 
gies, C,{z) = ^i^scdm/'^Acdm? for cach cluster in the sample. We then multiply Kaastra et al's (2004) 
values for hh by a/^ (since the observed X-ray flux and angular size are fixed) and multiply linear 
distances by C,^^ . Values of ^, as well as the cooling radius, are given in table |2l [The conversion 
from riH to rig is given by equation (|2T]) .] 

We adjust a single parameter, rsource? to fit to the observations. The optimal values for rsource 
are given in table [3l The temperature and density profiles in the model solutions are plotted in 
figures \T\ and |2l We note that the central peak in the Perseus temperature data is due to the hard 
power-law spectrum of the central active galaxy NGC 1275 (E. Churazov, private communication). 
The values of M and Lcr, as well as fluid quantities at ri = 0.2 kpc and router are given in table[3l The 
radial profiles of the rms turbulent velocity and cosmic-ray pressure are presented in appendix lAl 

Overall, the model profiles match the observations quite well, which suggests that convection 
is an important process in cluster cores. The model, however, substantially underestimates the ob- 
served density in Sersic 159-03 at r < 50 kpc. This discrepancy may be explained by a recent study 
by Werner et al (2007). These authors found that Sersic 159-03 has the largest soft x-ray excess of 
all clusters observed by XMM-Newton and argued that the observed excess is best explained by 
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Table 2: Cluster parameters 



Cluster 


^scdm 
^Acdm 


^ cool 

(kpc) 


Virgo 


1.0 


73 


Abell 


1.39 


61 


Sersic 159-03 


1.36 


128 


Abell 4059 


1.37 


86 


Hydra A 


1.36 


130 


Abell 496 


1.38 


89 


Abell 1795 


1.36 


130 


Perseus 


1.39 


128 



The quantity Jscdm is the angular-diameter distance to each cluster in the SCDM cosmology employed by Kaastra et al 
(2004) in which Q.= 1 andHf, = 50 km s^'Mpc^^ t/Acdm is angular-diameter distance to each cluster in the ACDM 
cosmology assumed in this paper, in which Q. — 0.7, Q.\ = 0.3 and Hq = 70 km s^'Mpc^'. For Virgo, t/scdm/^^Acdm is 
set equal to 1.0, since we use the same distance (16 Mpc) to Virgo as employed by Kaastra et al (2004). The cooling 
radii rcooi are taken from Kaastra et al (2004) butrescaled to ACDM. Kaastra et al's values of rcooi are the radii at which 
the radiative cooling time ?cooi is 15 Gyr in SCDM. At the rescaled values of rcooi. ?cooi is 15 Gyr-(JAcdm/<iscdm)^''^ 
in ACDM (because the cooling time scales like and He °^ d~^l^ for a fixed observed X-ray flux and angular size, 
where d is the angular-diameter distance). 

the presence of a substantial population of non-thermal electrons that is concentrated in the cluster 
core. When they modeled the observed emission as coming from a combination of thermal plasma 
and nonthermal electrons, they found that the x-ray emission between 0.3 and 10 keV from non- 
thermal electrons is a substantial fraction of the emission from the thermal plasma at large radii 
(~ 35 — 55% at r ~ 375 kpc) but only a small fraction of the emission at small radii (~ 1 — 7% at 
r < 50 kpc). Thus, the non-thermal contribution to the emission measure does not lead directly to 
a large change in the observationally inferred electron density at r < 50 kpc. On the other hand, 
Werner et al (2007) note that if there is a non-thermal proton population with significantly more 
pressure than the non-thermal electrons, the total cluster mass may be significantly underestimated. 
Moreover, since the non-thermal pressure inferred by Werner et al (2007) peaks strongly towards 
the cluster's center, the actual gravitational acceleration in the central 50-100 kpc may be much 
larger than in an NFW profile calculated neglecting non-thermal pressure. This is an issue for all 
the clusters that we consider, but especially for Sersic 159-03, since its especially large soft excess 
indicates a large non-thermal pressure fraction. We note that although the non-thermal emission 
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is less peaked than the thermal emission in the results of Werner et al (2007) (i.e., Pnon-thermai/^e 
decreases towards the center), the non-thermal pressure is more peaked than the thermal pressure 
(/'non-thermal/ P increases inwards). If we were to re-calculate our model solutions using a larger 
gravitational acceleration in the cluster core, the thermal plasma density would peak more sharply 
near the cluster center than in Figure [U Thus, the deviation between the model and observations 
of Sersic 159-03 may be due to a significant underestimate of the gravitational acceleration in this 
cluster resulting from its unusually large non-thermal pressure. 

Table 3: Physical quantities in the model solutions 



Cluster 


^ source 

(kpc) 




(keV) 


Pcrin) 

pin) 


l(v.(n))| 
cM) 


A/Bondi 

(Mgyr-i) 


(erg/s) 


^ outer 

(kpc) 


(keV) 


'^e, outer 

(cm-3) 




Virgo 


1.70 


0.115 


1.27 


0.268 


0.00207 


0.00203 


5.77 X 10^1 


49.3 


2.50 


2.85 X 10- 


-3 


Abell 262 


34.0 


0.231 


0.118 


0.153 


0.0388 


0.0232 


6.59 X 10^2 


93.6 


2.16 


1.22 X 10- 


-3 


Sersic 159-03 


46.0 


0.397 


0.133 


1.24 


0.365 


0.399 


1.13 X lO^"* 


165 


2.38 


1.35 X 10- 


-3 


Abell 4059 


24.0 


0.168 


0.727 


0.156 


0.0562 


0.0607 


1.72 X 10"*^ 


137 


3.89 


1.75 X 10- 


-3 


Hydra A 


35.0 


0.742 


0.340 


0.740 


0.257 


0.839 


2.38 X lO"*-* 


160 


3.28 


2.05 X 10- 


-3 


Abell 496 


5.00 


0.0437 


0.965 


0.701 


0.0201 


0.00650 


1.85 X 10^2 


96.2 


3.93 


2.53 X 10- 


-3 


AbeU 1795 


40.0 


0.671 


0.165 


0.700 


0.177 


0.363 


1.03 X 10^"^ 


184 


5.56 


2.19 X 10- 


-3 


Perseus 


16.0 


2.47 


0.487 


0.121 


0.0264 


0.343 


9.75 X 10^^ 


162 


5.27 


2.16 X 10- 


-3 



'"source is the size of the cosmic -ray acceleration region that leads to the best fit between the mixing-length model and the 
observations of Kaastra et al (2004). ng{ri), T{ri), p{ri), and Pa{ri) are the electron density, temperature, thermal 
pressure, and cosmic-ray pressure at the inner radius ri = 0.2 kpc. v,- and Cs are the radial velocity and adiabatic 
sound speed. Meondi is the Bondi accretion rate based on the plasma density and plasma temperature at r = ri , and 
La- = 0.005MBondic2 is the cosmic-ray luminosity of the central radio source. Touter and 'le. outer fhe observed 
temperature and electron density at the radius router of the outer boundary used in our shooting method. 

One of the quantities that is calculated as part of our solutions is the cosmic-ray luminosity Lcr 
of the central AGN. To further test the plausibility of our model, we compare the theoretically 
predicted values of Lcr from tableland the observationally inferred mechanical luminosity Lmech 
of the central AGN in the six clusters in our sample for which we were able to find published 
values. The mechanical luminosities are taken from Birzan et al (2004), and are calculated from 
observations of X-ray cavities, by assuming that an energy pV (where p is the surrounding pressure 
and V is the cavity volume) per cavity is released during a time equal to the buoyancy time scale. 
We plot Ljnech versus Lcr in figure[3l The error bars in this figure take into account projection effects 
on the estimate of the cavity volume as well as uncertainties in the ages of the cavities. These two 



- 16- 




r (kpc) r (kpc) r (kpc) r (kpc) 



Fig. 1 . — The solid lines give the electron density as a function of radius in our model solutions for the eight clusters 
in our sample. The data points are from the observations of Kaastra et al (2004). 

luminosities are strongly correlated over a range of ~ 100 — 1000 in luminosity. 



4. The energy budget of the intracluster medium 

In our model, radiative cooling is balanced by a combination of thermal conduction, convec- 
tive heating, and radial inflow due to the accretion onto the central AGN. To distinguish between 
these last two mechanisms, we separate the average radial velocity into two components, 

{Vr) = Vinflow + Vr,turb, (36) 

where 

- ^ 

is the inflow rate that arises in a laminar radial flow with constant mass accretion rate M. The 
term v^turb is an additional average radial velocity that is induced by the convection. [Its value is 
given by — (6p5v^) / (p), as in equation (|B5I ).1 With this definition in hand, we write the average of 
equation (HI), divided by (y — 1), as 

= (-^inflow + -f^conv + -f^visc + /^tc — R) ■ (38) 




Fig. 2. — The solid lines give the temperature as a function of radius in our model solutions for the eight clusters in 
our sample. The data points are from the observations of Kaastra et al (2004). 

Here, 

(^inflow) = - (^_\)^2 ^ {r\naoM) " {r\nflo^) (39) 

is the source term associated with Vinflow The term 

(^conv) = (-^^^T^-pV-^-^mflow^ (40) 

is the convective heating rate of the thermal plasma, excluding viscous dissipation. It includes the 
turbulent diffusion of heat as well as the turbulent pdV work done on the thermal plasma by cosmic 
rays. The average of the viscous dissipation term is set equal to 

0.420^3 

(//visc) = (41) 

where I = OAr is the mixing length, Mnns is the rms turbulent velocity defined in equation (IB71I) . 
and the constant 0.42 is taken from direct numerical simulations of compressible magnetohydro- 
dynamic turbulence (Haugen, Brandenburg, & Dobler 2004)0 The average of //tc is given by 
equation (|26|) . and the average of R is given by equation (fTTI) . 



'The constant 0.42 is obtained by taking the mixing length / to correspond to n/kp in the simulations of Haugen 
et al (2004), where is the wave number at which kE{k) peaks, and E{k) is the power spectrum of the turbulent 
velocity. 
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Fig. 3. — Comparison between the cosmic-ray luminosities in our model calculations La- against observationally 
inferred values of the mechanical luminosities of the central AGN in six of the clusters in our sample (Birzan et 
al 2004). The dotted line represents equaUty between these two quantities. 

In figure m we plot the averages of //inflow (dotted line), //conv (long-dashed line), Htc (short- 
dashed line), and R (solid line), integrated over volume from the inner radius of our model (ri = 
0.2 kpc) out to radius r. We find that (//vise) ^ (//cony) everywhere in each cluster, and so we omit 
//vise from the figures to keep the plots easier to reado 

We can divide the clusters into two groups, those with AGN-dominated heating and those 
with conduction-dominated heating. In Hydra A and Sersic 159-03, the heating within the cen- 
tral 100 kpc is dominated by convection driven by the central AGN. In all of the other clusters, 
the heating within the central 100 kpc is dominated by thermal conduction. The inability of con- 
duction to balance cooling in Hydra A and Sersic 159-03 was previously noted by Zakamska & 
Narayan (2003). These authors constructed density and temperature profiles for clusters assuming 
that radiative cooling is balanced by conductive heating, setting the thermal conductivity equal to a 
constant fc times the Spitzer thermal conductivity. For Hydra A and Sersic 159-03, they found that 
the values of fc that best fit the observations were 1.5 and 5.6, respectively, much larger than the 
theoretically expected value of fc ~ 0.1 —0.2. In contrast, for Abell 1795 the best-fit value of fc 
was 0.2. 



^We note that equation (|38] l is not exactly satisfied by our model solutions. In our model, we use the total-energy 
equation [equation ( IB3b 1 instead of the plasma energy [equation ©I. Although equation (01 is exactly satisfied when 
equations ( IB3l l. (O, and (|5]l are satisfied, our mixing-length approximation of the average of equation (|4| is not exactly 
satisfied when our mixing-length approximations to the averages of equations ( IB3b . OJ, and Q are satisfied. This 
discrepancy is noticeable at the largest radii in Hydra A and Sersic 159-03. 











Hydra A 












// 






//'-■, 








/■'' 
/-'' 1 








/' ' ' 
/ ' 
/ / ' 
' / 1 







10 100 





10 100 



r (kpc) 



r (kpc) 



r (kpc) 



r (kpc) 



Fig. 4. — The energy sources and sinks in the thermal plasma, integrated over volume from the inner radius r\ = 
0.2 kpc out to radius r. The solid line is the radiative losses (X-ray luminosity), the short dashed line is the heating 
from thermal conduction, the dotted line is power contributed by the inflow associated with the mass accretion rate, 
and the long-dashed line is the heating power due to the convective turbulence, which includes both the turbulent 
diffusion of heat and pdV work by cosmic-rays. 



What are the characteristics of a cluster that determine whether AGN feedback or thermal 
conduction is the dominant heat source within the central 100 kpc? Given that K5 T^/^, it seems 
clear that a lower average temperature Tavg makes thermal conduction less able to balance radiative 
cooling, leading in turn to a relatively greater role for AGN feedback. However, although Tgyg is 
important for determining the relative strength of conduction and AGN feedback, on its own the 
value of Tavg does not explain our results, since Virgo and Abell 262 have temperatures comparable 
to that of Sersic 159-03 and lower than that of Hydra A. An equally important factor appears to 
be the baryon density in the core. In particular, the clusters with AGN-dominated heating are 
significantly denser at a given radius than clusters with conduction-dominated heating that have 
similar average temperatures. For example, for radii between 20 and 50 kpc, the electron density 
in Sersic 159-03 is 2-3 times greater than in Virgo. Similarly, at r = 100 kpc the electron density 
in Hydra A is ~60% larger than in Abell 4059o even though Abell 4059 has a higher average 
temperature and larger virial mass. (Typically, at a fixed radius the density is larger in hotter. 



^This ratio is based on the model density profile due to the offset in the radius of the observed densities, but a 
similar conclusion is reached by interpolating between data points. 
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more massive clusters.) It thus appears that the clusters in which AGN feedback dominates most 
strongly over conduction are those in which the clusters' ongoing formation channels unusually 
large quantities of baryons towards the clusters' cores. 

We conclude this section with a few additional comments relating to figures |3] and |4l Birzan 
et al (2004) found that their observationally inferred values of L^ech for 16 clusters were corre- 
lated with the X-ray luminosity inside the cooling radius rcooh denoted Lx, supporting the idea that 
AGN feedback is at least part of the solution to the cooling-flow problem. However, the level of 
the mechanical luminosity in their study turns out to be a factor of 1 to 20 lower than the X-ray 
luminosity, which raises the question of whether the mechanical luminosity is sufficient to offset 
cooling in these clusters. Our model solutions and figure [3] show that the mechanical luminosity is 
indeed sufficient when thermal conduction is also accounted for, and our discussion above regard- 
ing AGN-dominated heating versus conduction-dominated heating offers an explanation for the 
large variations in the ratio Lx/Lj^ech- When heating is dominated by conduction, AGN feedback 
heating is only a small fraction of the power radiated from within the cooling radius, and Lx/L^ech 
is large. On the other hand, when AGN-driven convection dominates the heating, the AGN heating 
power is similar to the total power radiated from within rcooi, and Lx ~ ^ech- 

We note that figure |4] shows that in Hydra A, Sersic 159-03, and Virgo, conduction actually 
acts to cool the plasma over a limited range of r due to the local maximum in the temperature 
profile. Also, in none of the clusters is the total convective heating rate equal to L^. This is 
because the cosmic-ray luminosity is the power deposited into the cosmic-ray fluid, and only part 
of this is transferred to the thermal plasma through pdV work. Additional plasma heating arises 
from the redistribution (turbulent diffusion) of plasma thermal energy resulting from the convective 
motions. As can be seen from figure |4] and table [3l the total convective heating of the thermal 
plasma is typically on the order of one-third of L^r- 

5. Clusters with central cooling flows 

An important point to emerge from figure |4] is the appearance of central cooling flows in 
several clusters - Abell 262, Sersic 159-03, Hydra A, Abell 1795, and Perseus - with radii r^f 
that are typically a few kpc. Heating is unable to balance cooling at r < r^f in these clusters for 
several reasons: the AGN feedback heating is distributed over a large volume, thermal conduction 
becomes less efficient at small r due to the lower temperatures and the fact that Kj- T^/^, and 
the radiative losses per unit volume peak sharply at small r due to the large plasma densities. As 
a result, a cooling flow develops in which the energy lost to cooling is replenished by the inflow. 
Within this central region, we have the approximate relation ?cooi ~ ^inflow^ where ^inflow = ''/|vinflow| 
is the inflow time and ?cooi = l.SnkBT /R is the local cooling time. Because ^inflow ~ ^cooh the mass 
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accretion rate at r = rcf is approximately given by 



M~47rr^fp(rcf) 



.^cool('"cf) J ^cool('"cf) 



(42) 



where is the radius of the central cooling flow region, and M^f is the mass of plasma contained 
within the cooling flow region. Because we have no sources or sinks of plasma, M is independent 
of r in our model. 

How do these central cooling flows match onto adiabatic Bondi flow at smaller radii? In our 
model, as r decreases from rcf towards zero, p rises and T decreases until the Bondi accretion rate at 
the fixed radius ri = 0.2 kpc matches the cooling-flow mass accretion rate given by equation (|42l) . 
However, our forcing the flow to become adiabatic at ri is artificial, and leads to an unrealistic 
plasma profile near ri with an abrupt transition in the flow at ri . A better approach was adopted by 
Quataert & Narayan (2000). These authors investigated radial inflow with cooling in the absence 
of thermal conduction and cosmic rays using a numerical shooting method and solved all the way 
in to the sonic point, r = rgonicj at which Vr = — c^. They found a smooth transition from an outer 
cooling flow with ^inflow — ^cooi to an inner adiabatic Bondi flow with ^inflow ^ ^cooh provided that 
^sonic < nr, where 

ru^^= 0.05 kpc J [soo^j (43) 

is the radius within which gravity is dominated by the black hole and a is the circular velocity 
of the BCG, which was taken to be independent of r. They also found that equation (|42l) was an 
accurate estimate of the numerically calculated mass accretion rate in the absence of mass dropout. 
For rsonic < hr, Quataert & Narayan's solution satisfies Cj ~ o = constant at r > rtr. In the absence 
of mass dropout, the condition ?cooi = ^inflow leads to the relation p oc ^^3/2 within the cooling-flow 
part of their solution. The Bondi accretion rate Msondi [given by equation (l30l) 1 evaluated at a 
radius r within the cooling-flow part of their solution thus increases towards smaller r like r^^/^. 
At a sufficiently small value of r, which we call rad, the Bondi accretion rate equals the rate Mcf 
at which mass flows in through the cooling flow, and the flow makes a transition to an adiabatic 
Bondi flow. At r < rad, the ratio ?cooi/^inflow increases towards smaller r, and so the neglect of 
cooling at r < rad is self-consistent. We note that the Bondi accretion formula can be applied in the 
model of Quataert & Narayan (2000) at the outer boundary of the adiabatic flow region, rad, even 
if rad lies outside the region in which the black hole dominates the gravitational acceleration. This 
is because the Bondi accretion rate depends only on the specific entropy s of the plasma and Mbh, 
and s is constant for r < rad- 

It would be valuable to incorporate into our model an approach similar to that of Quataert 
& Narayan (2000), including cosmic rays, thermal conduction, and the possibility of convection. 
Although such a calculation is beyond the scope of this paper, we expect that in such an analysis 
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the mass accretion rate of the central accretion flow and the plasma parameters at rBondi still 
controlled by Mcf, as in Quataert & Narayan's (2000) work. Because the central accretion flow is in 
some sense slaved to the surrounding cooling flow, the AGN regulates the mass accretion rate pri- 
marily by controlling the properties of the central cooling flow, and in particular by controlling r^f. 
For example, if M rises above the equilibrium value, the AGN-feedback heating rises. This then 
reduces rcf, because one has to go to smaller r in order for p to rise enough that cooling exceeds 
the convective heating rate. The reduction in r^f reduces M, as can be seen from equation (|42l) . 
which then causes M to drop back down to its equilibrium level. 

As mentioned above, the existence of a smooth transition from a cooling flow to an inner 
adiabatic Bondi flow requires that rsonic < hv- Otherwise, as described by Quataert & Narayan 
(2000), the ratio ?cooi/^inflow decreases inwards in the supersonic region at rtr < r < rsonic, and the 
plasma cools rapidly to very low temperature. In this case, the cooling plasma could still end up 
fueling the central black hole, but it would do so through some process other than the one we have 
assumed in our model, e.g., by forming stars whose winds then feed the black hole or through 
infalling cold gas [see, e.g., Pizzolato & Soker (2005) and Soker (2006)]. 

Under what conditions is rsonic > '"tr? One factor that can cause rsonic to exceed rtr is a small 
central black hole mass, since a smaller Mbh reduces r^. In addition, as illustrated in Quataert & 
Narayan's approximate analytic results, if the black hole's contribution to gravity were hypotheti- 
cally ignored, rsonic would increase with increasing M. Thus, a sufficiently large M can also cause 
''sonic to exceed r^. A large M results from either a large Lcr or a small accretion efficiency r\. 
As discussed in the previous section, is approximately determined by the baryon density and 
temperature at the cooling radius - a higher density and/or lower temperature at rcooi means that 
thermal conduction can offset less of the radiative cooling within the cooling radius, which in turn 
leads to a larger L^. Thus, to summarize, smaller values of Mbh, ^1, or r(rcooi) and/or larger values 
of p(rcooi) can cause rsonic to exceed rtr, preventing a smooth transition from a cooling flow to an 
inner Bondi flow, and causing the cooling of intracluster plasma to low temperatures at r > rtr. 

We note that if the accretion efficiency r|cool that arises when plasma cools to low tempera- 
ture outside rtr is much smaller than the accretion efficiency r\ associated with Bondi accretion, 
then a flow that cools rapidly outside rtr will need a much higher M (and larger rcf) in order for 
AGN feedback to provide the heating needed to offset cooling within the cooling radius. A much 
larger M, in conjunction with plasma cooling to low temperatures outside rtr, would imply a much 
larger star formation rate within the BCG. Thus, a small Mbh or large Lcr could lead to the condi- 
tion rsonic > fti and cause star formation at rates that significantly exceed the Bondi accretion rates 
listed in table |3l 

Returning to our model calculations, we list in table [3] the values of | (v, ) l/c^ at r = ri for the 
clusters in our sample. The value of | (v,) l/c^- reaches its maximum at r = ri in our solutions, and 
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thus our solutions satisfy \{vr)\/cs < 1 at all radii. Our calculations are thus at least marginally 
consistent with our assumptions of hydrostatic equilibrium and Bondi accretion. We also note that 
the the sound crossing time is shorter than cooling time t^ooi at all radii in our model solutions. 
However, our model imposes an abrupt and artificial transition in the flow at r = 0.2 kpc, which 
causes our solution near r = 0.2 kpc to be inaccurate. In addition, there is significant uncertainty 
in the values of Mbh and r|. It is thus possible that some of the clusters reach a sonic transition 
outside rtr. Further investigation of this issue is needed. 



6. Summary 

There is a growing consensus that AGN feedback holds the key to solving the cooling-flow 
and overcooling problems for clusters of galaxies. However, the way in which an AGN's power is 
delivered to the diffuse intracluster plasma is still not well understood. In this paper, we suggest 
that an AGN's mechanical luminosity heats the intracluster plasma by accelerating cosmic rays 
that cause the intracluster medium to become convectively unstable. We explore this idea by 
developing a steady-state, mixing-length-theory model. By adjusting a single parameter in the 
model (the size of the cosmic-ray acceleration region, rsource)> we obtain a good match to the 
observed density and temperature profiles in seven out of the eight clusters in our sample. Our 
model underestimates the density in the eighth cluster, Sersic 159-03, within the central ~ 50 kpc. 
We suggest that this discrepancy may result from the fact that the parameters in our NFW mass 
model are determined neglecting the cosmic-ray pressure. At the same time, Sersic 159-03 has the 
largest soft x-ray excess of any cluster observed by XMM, and likely contains a large population 
of non-thermal particles concentrated in the cluster core. (Werner 2007) If the mass model were 
recalculated taking the non-thermal pressure into account, the gravitational acceleration would be 
larger, especially in the cluster core, which would increase the plasma density at r < 50 kpc in 
our model calculations and possibly bring the model into agreement with the observations. We 
also find that the cosmic-ray luminosities of the AGN in our sample are strongly correlated with 
the observationally inferred mechanical luminosities of these AGN. Our results suggest that AGN- 
driven convection is an important process in cluster cores. 

In our model solutions, the radiative cooling rate is much more peaked about r = than is 
the rate of convective heating. As a result, a compact central cooling flow arises in our model 
calculations for several of the clusters in our sample. The radii, r^, of the cooling flows are 
typically a few kpc. The mass accretion rate onto the central AGN in these clusters is roughly 
the plasma mass at r < rgf divided by the cooling time at r^f. We suggest that the AGN regulates 
the mass accretion rate in these clusters by controlling rgf: if the AGN power rises above the 
equilibrium level, the size of the central cooling flow decreases, the mass accretion rate drops, and 
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Fig. 5. — The rms turbulent velocity as a function of radius in the model solutions. 

the AGN power then drops back down to the equilibrium level. 
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A. The profiles of the cosmic-ray pressure and turbulent velocity 

The profiles of the rms turbulent velocity Wnns [defined in equation (IB71I) 1 and the cosmic- 
ray pressure (as a fraction of the thermal pressure) are plotted in figures [5] and [6l Although it is 
difficult to see in the cases of Sersic 159-03, Hydra A, Abell 496, and Abell 1795, figure[6] shows 
that {d/dr){pcx/ p) > at small r. [This can also be seen by comparing the figures with the values 
of Pcx{ri) / p{r\) listed in table[3l] This is not because dpc^/dr > (in fact, dpcv/dr < at all r 
for each cluster), but instead because the thermal pressure decreases with radius more rapidly than 
the cosmic-ray pressure. We note that there was an error in one of the plotting subroutines used for 
paper II, which resulted in the velocities plotted in figure 4 of paper II being too large by a factor 
of 3. 
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Fig. 6. — The ratio of cosmic-ray pressure to thermal pressure in the model solutions. 

B. The two-fluid mixing-length theory and numerical method 

In this appendix we describe the two-fluid mixing length theory that we use to model the 
convective intracluster medium. We take all fluid quantities to be the sum of an average value and 
a turbulent fluctuation, so that 

p = (p)+5p, (Bl) 

v=(v)+5v, (B2) 

etc. We take all average quantities to depend only on the radial coordinate r, and we set (v) = {vr)f. 
We wish to solve for four quantities: (p) , (v^) , {T), and (pcr) • To do so we take the averages of four 
equations: equations ©J ©, and ©, as well as the total-energy equation. The latter is obtained 
by taking the dot product of equation (|3]) with v and adding the resulting equation to the sum of 
equations dH) and (|5]), which yields: 

d ( pv^ ^ p Per 
dt\ 2 y-1 Y„- 

^ / pVV^ ^ Jvp YcrVPcr ^ D ■ Vpc 

+v- ^ + pv4>+-?^ + -^^^^+rvi,c-K-vr ^ 

V 2 y-1 Ycr-l Ycr-l 

— P-^^^ ^+£^source, (B3) 
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where Fyisc is the viscous energy flux, and where we have made use of the relation //vise — (V • 
nyisc) ■ V = —V • FviscH We assume that the viscous energy flux is much less than the advective 
energy flux and thus drop Fyisc- We also set d^/dt = 0. 

The average of equation Q can be written 



1 

(P) 



Q + 



M 



(B4) 



where 



e=(5p6v,), (B5) 
and the mass accretion rate M = — 47rr^(pv, ) is a constant. The average of equation ([3]) yields 



ar 



-(p)-r' 

ar 



where 



Pm = P + Pc 



(B6) 



(B7) 



is the total pressure. In writing equation (|B6I) . we have taken the convection to be subsonic, so 
that the Reynolds stress can be neglected. We have also dropped the viscous stress, which is 
unimportant in the averaged equation. The average of equation ([5]) can be written 

d^Pcr) 



D„- 



(1-Y) d ,2 



Jr2 



dr 



(rV)+(l-Ycr) (W+(£source)) 



djpci) 
dr 



1 J / \ 
r^ dr 



+ 



[r (v,)) + (v,) 



where 



and 



dr 



F = 



djpcv) 
dr 



Y-1 



(B9) 
(BIO) 



W = {dpV-dv). 

In writing equation (|B8I) . we have again made use of the fact that the convection is subsonic, 
which implies that the total-pressure fluctuation is very small. As a result, we can set dp ~ —dpci, 
(5v^5j9cr) = —{dvr^p), and (Spcr^ ■ 5v) = — (5pV ■ 5v). The average of equation (IB3I ) yields 

d^T) 



Kt- 



dr^ 



d^ 

Q— + 



yp 



dr {^—\)r^dr 



Y-1 



(i4> dp a 
^ dr dr 



^ 1 J . 2^\ 

r^ dr 



^The fact that the viscous terms can be written as a total divergence reflects the fact that viscosity is neither a 
source of energy nor a sink of energy, but instead merely converts bulk-flow energy into thermal energy. Thus, when 
equation ( |B3t is integrated over volume, Gauss's law can be used to express the viscous terms as a surface integral, 
which vanishes if the boundary of the integration lies outside the plasma. 
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dr 



(BID 



where we have dropped the (pvv^/2) term since it is much smaller than the other terms for sub- 
sonic convection. In equation (IBl II) . is the average (isotropic) thermal conductivity given by 
equations (|24l) and (l25l) . where T is set equal to (T) in equation (|24)) . 

To solve equations dHU), (Ell), dBS]), and (IBllI) for (p), (v^), (r), and (pcr), we need to 
express the quantities Q, F, and W in terms of these average fluid quantities, thereby closing the 
equations. We accomplish this by using a two-fluid mixing length theory. Our approach is to first 
estimate Q, F, and W using a local mixing length theory. In the local theory, the properties of 
the turbulence at some radius are determined only by the average fluid properties and gradients at 
that radius. We then use this local theory as the basis for a non-local mixing length theory, as in 
paper II. In the nonlocal theory, the properties of the turbulence at some radius r are determined 
by a weighted average of the turbulence properties in the local mixing length theory over a range 
of radii. Our mixing-length theory differs from stellar mixing-length theory (Cox & Giuli 1968) 
in two important ways: we include a cosmic-ray fluid, and we take into account the fact that the 
diffusion of heat and cosmic rays occurs almost entirely along magnetic field lines. 

We derive quantities in the local mixing length theory as follows. We take the convective 
turbulence to have a correlation length Z, also called the mixing length, where 

/ = ar, (B12) 

a is a constant, and r is the distance from the center of the cluster. Fluid parcels in convective 
regions are taken to rise or sink a distance / before breaking up and mixing into the surrounding 
plasma. We take / to be much smaller than the pressure scale height, so that a is treated as ^ 1 . 
(However, as in standard mixing-length theory, after the mixing-length-theory equations are de- 
rived, we relax the requirement that a <^ 1, and set a = 0.4 when applying the model to clusters in 
section|3l) We treat the fluctuations as small quantities, and take a/ (5p)^^ /(P)' a/ (Sr)^^ / (T), 

<^a/(5p^^^ / {pcr), and <^a/|5vP^ /cg to be ~ o(a) (meaning of order a), where Cg is the sound 
speed. We then expand the equations in powers of a, and keep only the lowest-order non-vanishing 
terms in this expansion. The quantities Q, F, and W involve products of fluctuating quantities, and 
are thus ~ o(a^). We take R, (^source), Kj, K||, Dcr, £>||, and //vise to be ~ o(a^), so that, e.g., 
radiative cooling, conduction, and turbulent heating are of the same order in a in equations (|B8I) 
and (IBl II) . (Since the fluctuations are small, we can write, e.g., that (p) ~ kB{p){T) //uniH, where 
/J is the mean molecular weight.) There are two contributions to the average velocity (v). One 
is driven by the turbulent fluctuations, and, as will be seen below, is of order a^. The second 
arises from the net inflow of mass towards the center of the cluster. In our model, the mass ac- 
cretion rate is set by the Bondi accretion rate calculated from the plasma parameters at the inner 
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radius ri = 0.2 kpc, as described in section [2l We treat this second contribution to (v) as also of 
order a^. 

We now proceed to estimate typical values for dT, 6p, and dpa to first order in a - i.e., 
ignoring terms of order a^. Because we take the fluid displacement and the correlation length to 
be comparable, we need to use a Lagrangian approach to formally integrate the equations. We take 
the initial position of a fluid element at time ? = to be denoted ro, and its position at time t to be 



(B13) 



where ^ is the displacement of the fluid element. We use the shorthand notation that p{t), T{t), 
p{t), Pcx{t), and ^(?) are the density, temperature, pressure, cosmic-ray pressure, and displacement 
at time t of the fluid element that started at position tq dXt = 0. The velocity at time t of the fluid 
element that starts at ro at ? = is given by 



dt 



We then have that 



a 



d_ 



dt 



+ V-V, 



(B14) 



(B15) 



where the spatial derivatives on the right-hand side are with respect to r, not tq. The Jacobian 
matrix 7 for the transformation from ro to r is given by the equation 



Ji 



droj' 



(B16) 



The determinant of this matrix, denoted J, satisfies the equation 

= V-v, 



|-ln7 
dt 



(B17) 



'■0 



where the spatial derivatives on the right-hand side are again with respect to r, not Using 
equations (IB 151) and (|B17I) . we rewrite equation ^ as 



ainp 



ro 



din J 



(B18) 



^Equation IB 1 71 can be shown as follows. We define M,,- to be the determinant of the 2x2 matrix obtained by 
deleting the /''^ row and column of the matrix 7. We can then write that dJ /dJij = (— l)'+^M,j. The inverse of 
J_, which we denote B, satisfies BijJjk = 8,i, and is given by Bij = (— l)'+'My,/y. We can thus write dJ /dt\^^ = 
{dJ/dJij){dJij/dt)\ = JBji {dJij/dt) | . From the chain rule, (8ro,/3ry) (3r,73roit) — 5,^. Thus, Bij = droi/drj, and 
(l/J) dJ/dt\,.^ = idroj/dn)[id/dt)\,.^ (dn/droj)] = {droj/dn){dvi/droj) = dv./dr,. 
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We integrate equation (IB 181) in time from t = 0to t = At holding ro fixed (as in the time integrals 
below), where 

At = — (B19) 

ML 

is the "mixing time," and ul is the rms radial velocity in the local mixing length theory. We set 
^(0) = and p(0) = (p), where (p) is shorthand notation for the average density at tq. We then 
obtain 

p(AO 1 

where J{At) is the Jacobian att = At evaluated for the initial position tq. If we start at ? = with 
a fluid element of infinitesimal volume d^ro centered at ro, then at time At its volume is d^r = 
J{At)d^rQ. Thus, equation (IB20|) is a statement of mass conservation. The Lagrangian density 
perturbation at time At is 

ApLag = P(A0 -P(0) = p{At) - (p). (B21) 
Combining equations (IB20I) and (IB21I) . we write 

To solve for the pressure fluctuation, we write equation @ in the form 

dp 



dlnJ 

IP 



+ (Y-l)5//tc, (B23) 



dt 

where 

i/tc = V-(K||S^-Vr) (B24) 

is the rate of heating due to thermal conduction, and 5//tc is the deviation of H^c from its average 
value. In writing equation (|B23I) . we have dropped terms of order a^. [5//tc ~ 0(a) since hT ~ 
0{a), K|| ~ 0{a?-), and W^hT ~ hT /f- ~ o{a^^).] We integrate equation (IB23I) from ? = to 
t — At, setting ^(0) = and p{0) = {p), where (p) is the average density at ro. To first order in a, 
we can replace yp d\nJ/dt\^^ with y{p) d\nJ/dt\^.^. We thus obtain 

A;?iag = -Y(;?)ln7(A0 + (Y-l) ^ dH,,{t)dt, (B25) 

Jo 

where 

ApLag = p{At)-p{Q) = p{At) - {p), (B26) 

is the Lagrangian pressure perturbation at time At. Because the density fluctuations are small, 
the value of J {At) is very close to 1. Writing J {At) = l-\-x, the quantity x is of order a. Thus, 
\nJ{At) =x+0{a?) = \-[J{At)]-^ + o{a?) and 

\nJ{At) = -^^ + 0{a}). (B27) 
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To first order in a, we can thus rewrite equation (IB25I) as 

Mag = I^^l^ + (Y- 1) l%H,,(t)dt. (B28) 

We have been analyzing a fluid element that starts off at ? = as an "average" fluid element, 
and we thus take 5r = and 5//tc = at ? = 0. When this fluid element is displaced radially 
outwards a distance /, it remains magnetically connected to the same set of fluid elements to which 
it was initially connected (at least until it is mixed into the surrounding fluid, at which point it is 
assumed that the magnetic field in the fluid parcel is randomized). This is depicted schematically 
in figure |7l If its temperature remains unchanged as it moves [i.e., if STLag = T{t) — T{0) = 0], 
then the effects of thermal conduction are unchanged and 6//tc will remain zero. However, if 
its temperature decreases, then more heat will be conducted into the fluid element, and 5//tc will 
increase. Thermal conduction will thus act to restore the temperature to its initial value (i.e. to 
keep 6rLag = 0), with 5//tc °^ — STLag for small values of dT^^g. To estimate 5//tc, we note that the 
mixing length / is comparable to the correlation lengths of both the convective turbulence and the 
temperature fluctuations. Thus, |V5rLag| ~ STLag//, and we have the order-of-magnitude relation 

6//tc (B29) 



We estimate that 



, f^^ 0.3(Y-l)K||ArLagA/' 

(Y-1)/ dH,,dt = ^-^ y , (B30) 







where ATLag = 5rLag(Af) is the Lagrangian temperature perturbation at time At, and the numerical 
factor of 0.3 is chosen somewhat arbitrarily to reflect (1) our expectation that the length scale of 
the temperature fluctuations is somewhat larger than Z, which is just the radial component of the 
displacement vector, not the full modulus of ^, and (2) the fact that 6//tc increases from zero to its 
maximum value as t ranges from to At, so the typical value of 6//tc is less than its value att = At. 
Finally, we obtain the relation 

Y(p)ApLag 0.3(Y-l)K||ArLagA? 
Aplag = J2 • (B31) 



To first order in a. 



ATLag _ A/?Lag ApLag ^^^2) 



{T) ip) (P) ' 
where (T) is the average temperature at tq. Equations (IB31I) and (IB32I) combine to give 

jp) f y+aj 
(P) \l+ai 



ApLag = ApLag^ ( ) ' (^33) 



-31 - 




Fig. 7. — Schematic diagram of a rising fluid parcel. The solid line is a magnetic field line passing through the 
parcel's initial location. The dashed line is an idealization of how the field line changes as a result of the fluid parcel's 
displacement. 



where 

Q-3(y-i)K||(r) 

ai = — —f (B34) 

is roughly the ratio of the mixing time At to the time for heat to diffuse a distance / along the 
magnetic field. When ai <^ I the thermal plasma expands adiabatically, and when ai ^ 1 the 
thermal plasma expands isothermally. 

To solve for the cosmic-ray pressure fluctuation, we write equation ([5]) in the form 
dp, 



dt 

where 



ain7 

YcrPcr 



+ 5//diff + (Ycr - 1 ) 5£source , (B35) 



Him = y-{D{bb-Vp,,), (B36) 

6i/diff is the deviation of //jiff from its average value, ^source is the cosmic-ray energy per unit 
volume generated by the central radio source, and S^somce is the deviation of ^source from its 
average value. As discussed in section 12.71 S^somce can be significantly larger than (^source) if 
the cosmic rays are accelerated in a small fraction of the volume. We thus treat S^source as O (a) 
and (^source) as C'(a^). We integrate equation (IB35I) from t = to t = At, setting ^(0) = and 
Pct{0) = (Pcr), where (pcr) is the average cosmic-ray pressure at tq. To first order in a, we can 
replace JcrPcr ^lny/3f|^y with Ycr(Pcr) ^lny/3r|^^. Using equation (IB27I) . we obtain 

ApcrXag = ^^^^^^^^ + 1^' [5//diff + (Ycr - l)5£source] dt, (B37) 
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where 

Apcr,Lag=i^cr(AO-(;?cr), (B38) 

is the Lagrangian cosmic -ray pressure perturbation at time A?. 

We treat parallel cosmic -ray diffusion in the same way as parallel thermal conduction and 
make the estimate 

rAf 0.3D||Apcr,LagA? 



/„ 



dHiiffdt = " ; ' ^ ■ (B39) 

We assume that S^source is typically positive in outwardly moving fluid elements and negative 
in inwardly moving fluid elements. Here, we are focusing on a a fluid element that is moving 
outwards (the discussion can be repeated with little alteration for inwardly moving fluid parcels), 
and thus we treat S^source as positive. We then make the estimate that 

/•At 

/ 8E,ouicedt = 8E,^,At, (B40) 
Jo 

where dE^^s is the rms value of S^source^ which is determined from equation (|35l ). Substituting 
(IB39I) . and (IB40I) into equation (IB37I) . we obtain 

A;.cr,Lag = ^^^YT^^ + (Ycr " ^nnsX. (B41) 

where 

'uL 0.3D 



is the effective time during which cosmic rays can accumulate in the fluid parcel as a result of the 
cosmic-ray source term, (x is roughly the shorter of the mixing time l/ui and the diffusion time 
.) The quantity a2 is given by 

0.3D|| 

ai = -r^ (B43) 
/Ml 

and is approximately the ratio of At to the time for cosmic rays to diffuse a distance / along the 
magnetic field. When a2 -C 1, the cosmic rays expand adiabatically if dEi^s = 0. When a2 ^ 1 the 
cosmic ray pressure in the fluid element remains constant as the element is displaced if S^rms = 0, 
as in the linear Parker instability in the large-D|| limit (Parker 1966, Shu 1974, Ryu et al 2003). 

Since it is assumed that the convection is subsonic, the total pressure in the fluid element 
remains approximately the same as the average total pressure in the fluid element's surroundings. 
We take 

r-m) = l: (B44) 

and thus to first order in a 

ApLag + Apcr,Lag = I ^ (i^tot) • (B45) 
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Adding equations (|B33I) and (|B41I) and using equation (IB45I) . we obtain 

d 



/^(;?tot) = CeffApLag + (Ycr - l)5£rmsX, (B46) 



where 

'l+ai \ {P) , ycr(Pcr) 



Ceff 



+ ■ 



1+aiJ (p) (l+a2)(p) 
is an effective sound speed for the medium. 



(B47) 



The fluctuating quantities appearing in equations (|B5I) . (IB9I ). and (|B10I) are Eulerian fluctu- 
ations, in that they involve the difference between some quantity and the average of that quantity 
at the same location. To first order in a, we can write the Eulerian density fluctuation of our 
outwardly displaced fluid element at time At, denoted Ap, as 

dip) 

Ap = ApLag-/-^. (B48) 
Combining equations (IB46I) and (IB48I ). we find that 

Ap^lf - ^ V - (B49) 



\c;ff dr dr / c;^ 

The fluid is convectively stable if an outwardly displaced parcel is heavier than its surroundings 
(i.e., if Ap > 0) for any value of ui. We note that as ml increases, ai and a2 decrease, Cg^^ increases, 
and X decreases. Since d{pioi)/dr < 0, it follows that 

^ Ap>0. (B50) 



dui^ 



Thus, if Ap > as ui 0, then Ap > for any {ui is by definition non-negative), and the 
medium is convectively stable. On the other hand, if Ap < as ui 0, then the medium is 
convectively unstable. The necessary and sufficient condition for convective stability is thus that 
Ap be positive in the limit ui 0. As ul 0, we have that ai — oo, a2 ^ °°, x ^ l^/ (0.3D||), and 
'-eff ~^ (p) I (P)- constant mean molecular weight ^u, this then leads to the stability criterion 

^ J7 + 03^ > ^' ^^^'^ 

where n = p/{ijm{j) is the number density of thermal particles. If one sets S^rms to zero, then 
equation (IBS II) reduces to the stability criterion derived by Chandran (2005) and Chandran & 
Dennis (2006). Here, we have kept the fluctuations in ^source, which act to destabilize the medium 
to convection, since localized excesses in the cosmic -ray pressure lead to pockets of buoyant, 
lower-density fluid. 
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If the convective stability criterion is satisfied at some radius, we set the local convective 
velocity ml to zero at that radius. Otherwise the fluid is convectively unstable, and we estimate the 
value of Ml by solving the polynomial equation 



(P)"l 



/Ap 



16 dr 



(B52) 



Equation (|B52I) states that the mean radial kinetic energy of the fluid element is the mixing length 
times the buoyancy force on the fully displaced parcel times the numerical factor of 1/16 that is 
commonly used in one-fluid mixing length theory (Cox & Giuli 1968). Once ml is found, we 
determine ApLag and ApLag using equations (IB46I) and (|B33I) . respectively. The Eulerian pressure 
perturbation at time At, denoted Ap, is then given by the equation 

Ap) 



Ap = ApLag-^ 



dr 



We then estimate the quantity F in equation (IB9I) to be 



-avg 



ui^Ap 



Y-1 



(B53) 



(B54) 



Here, as below, the L subscript is used to denote the estimate obtained using local mixing length 
theory. We set 

cavg = 1/2 (B55) 

to match standard treatments of local one-fluid mixing length theory (Cox & Giuli 1968). We 
estimate the quantity Q in equation (IB5I) to be 



Ql = CavgWLAp, 



(B56) 



with Ap determined from equation ( IB49I) . We estimate the quantity W in equation (IB 101) by noting 



/■At 

that / V-vdt 
Jo 

J_ / ApLa; 

At V (P) 



p{At) 
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(P) 

MApLag 



In 



1 + 



ApLag 
(P) 



^PLag 
" (P) 



. Thus, the typical value of V ■ 6v is 



. We thus set 



Wl 



MiApApLag 

l{9) 



(B57) 



Having estimated Q, F, and W using local mixing length theory, we now use these estimates 
as the basis for a nonlocal theory. In the study of Ulrich (1976), the nonlocal heat flux is given by 
a weighted spatial average of the heat flux obtained from local mixing length theory. We adopt the 
same approach and set 

/oo 
(B58) 
-oo 



-35- 



and 



/oo 
dziWi^{zi)^w{z-zi), (B59) 

/oo 
dziQi.{zi)y^Q{z-zi), (B60) 
-oo 



where 

z = lnf— y (B61) 

Tref is an unimportant constant, and the NL subscripts denote values in our nonlocal theory. Dif- 
ferent forms for the kernel function v^f were considered by Ulrich (1976). Here, we adopt the 
following values: 



¥eW=¥FW = <i Q if;c<0 ' ^^^^^ 



and 

ig-x/a^ if ;c > 



Equations (IB58I) through (IB63I) are equivalent to the differential equations 

ar— ^+FNL = ii, (B64) 
ar 

ar%5^ + eNL = eL, (B65) 
ar 

and 

aw'--P + WNL = WL. (B66) 
dr 

For r > Tconv, where rconv is the largest radius at which the fluid is locally unstable to convection, 
Fl = and Fnl °^ r^^/". When F and W are set equal to Fnl and V^nl in equations (|B8I) and (|B11|) . 
the terms containing Fnl are oc fQj. ^ y rconv To obtain the same scaling for the terms 

containing V^nl, the value of ttiy is determined from the equation 

a^''=a^^ + l. (B67) 

Equations (IB62I) and (|B63I) represent a one-sided average, in the sense that the nonlocal quantities 
^NL, 2nl, and Wnl depend only on the values of Fl, Ql, and Wl at smaller radii. A more sophis- 
ticated nonlocal theory could be developed along different lines (see e.g. Travis & Matsushima 
1973, Ulrich 1976, Xiong 1991, Grossman, Narayan, & Amett 1993), but is beyond the scope of 
this paper. 

The final equations for our mixing-length model are then equations (|B4I) . (|B6I) . (|B8I) . (IBl II) . 
(IB64I ). (IB65I ), and (IB66I ). which form a system of two second-order ordinary differential equations 
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(ODEs), four first-order ODEs, and one algebraic equation [equation (IB4I) 1 for the seven variables 
(p), (v^) (r), {pci), F^i^, 2nLj and W'nl- Eight boundary conditions are required to specify a so- 
lution. Two boundary conditions are obtained by requiring that the model density and temperature 
match the observed values pouter and Touter at the outer radius router- For seven of the eight clusters 
in our sample (all except Virgo), we choose router to be the center of the first radial bin outside the 
cooling radius rcooi- Values of rcooi for each cluster are given in tabled For Virgo, we take router to 
be the outermost data point, which lies inside of rcooi- Since we do not solve all the way in to the 
sonic point, we are forced to pick inner boundary conditions (at ri =0.2 kpc) in a somewhat arbi- 
trary way. We take d{T)/dr, d{pcY)/dr, Fnl, V7nl, and Qnl to vanish at r = r\. These "no-flux" 
boundary conditions set the diffusive and turbulent energy fluxes to zero at the inner boundary. 
Although this choice is undoubtedly inaccurate, we expect that it has only a small effect on our 
solution for the structure of the intracluster medium at r ^ ri . The eighth boundary condition is 
obtained by assuming that (pcr) ^ as r — >^ oo. This condition is translated into a condition on (pcr) 
at router as follows. The value of router is chosen to be significantly greater than rsource and much 
greater than Dq/vci, so that for r > router, ^'source is negligible and Dcr ^ v^rQ In addition, router is 
taken to lie outside rconv, the largest radius at which the intracluster medium is locally convectively 
unstable, so that Fnl = i^outer ('"/'"outer) "^^^ and H^nl = Wouter('"/'"outer)"^"^/'^ for r > router, where 
Fouter and Wouter are the values of Fnl and V7nl at r = router- We also take (v^) to be negligible 
for r > router- This latter assumption is reasonable, since Qnl = 2outer ( '"/'"outer) ^^'''^ for r > rconv, 
where Qouter is the value of 2nl at router- The resulting value of (v^) is significantly less than Vd 
for r > rconv in the numerical solutions we present in section |3l and thus (v^) plays only a small 
role in equation (IB8I ) at r > rconv- Solving equation (IB8I ) and requiring that (;?cr) ^ as r ^ oo, 
we find that for r > router 

d{Pcr) X 2(pcr) .-p.... 

dr Vdr^+v« r 

where 

X = (2a- 1)(Y- l)Fouterrler + «(ycr- l)Wouter'"iie/". (B69) 

Equation (IB68I) applied at r = router provides the eighth boundary condition. We then solve our 
system of equations using a shooting method. We guess the values of (p), (T), and (;?cr) at r = ri 
and then integrate the equations from r = ri to r = router- We then update our three guesses using 
Newton's method until the three boundary conditions at router are met. 

To compare to future observations and to analyze the turbulent diffusion of metals in the ICM 
[see, e.g., Rebusco et al (2005)], it is of interest to calculate the rms turbulent velocity. We define 



*The case Vd = requires a different approach and is not treated in this paper. 
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a nonlocal turbulent velocity, mnl, through the equation 

CLr — |-MNL = ML- (B70) 

ar 

Since ui (and thus mnl) is an estimate of the radial component of the velocity of a convective fluid 
element, the full rms turbulent velocity is roughly 

Mnns = \/3MNL, (B71) 

which is the quantity plotted in figure [51 
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